




*-------------------------------------------------------------------------------
* Start with crime data
*-------------------------------------------------------------------------------


*Start with 2008
use "offenses_known_monthly_2010.dta", clear

*Append other years
forval i = 2011/2015 {
append using "offenses_known_monthly_`i'.dta"
}

*Replace fips codes for NYC
replace fips_state_county_code = "1" if fips_state_county_code == "36005" ///
 | fips_state_county_code == "36061" | fips_state_county_code == "36047" | ///
fips_state_county_code == "36081" | fips_state_county_code == "36085" 

*Replace fips code for Richmond
replace fips_state_county_code = "51760" if fips_state_county_code == "51159"


*Prep variables
destring fips_state_county_code, gen(temp)
tostring temp, gen(Id2)
drop temp
drop month
gen month = date(date, "YMD", 2030)
replace month = mofd(month)
format month %tm



***
*merge in removals data
***

drop state
merge m:1 Id2 month using removalcounts, keepusing(firstmonth state county policymonth monthlysanctuary)

**drop pre-2010 and post-2015 data
drop if month < m(2010m1) | month > m(2015m12)

**We've now matched all observations in the removals data can drop excess crime data observations
keep if _m == 3
drop _m

**Now remove observations before first S. Comm. deportation
keep if month >= firstmonth



**
*Restrict to counties with 90% of crimes accounted for by never-missing agencies.
**
    
*mark those notes as having any missing observations
bysort ori: egen total_missing = sum(number_of_months_missing)

*now take a look at ones with unusual zeros
*Look at places with zeros in months; we need to drop places that don't report monthly
bysort ori: egen minori = min(actual_index_total)

*means by month
bysort ori: egen meanori =mean(actual_index_total)

*browse
sort ori month
br ori month actual_index_total month_mi if meanori > 50 & minori == 0 & total_missing == 0

**assume these are missing values if mean is 50 and there is ever a zero
replace total_missing = . if meanori > 50 & minori == 0 

*first total including all agencies
bysort Id2: egen countytotal = sum(actual_index_total)

*now total from nonmissing agencies
bysort Id2: egen temp = sum(actual_index_total) if total_missing == 0
bysort Id2: egen countynonmissing = mean(temp)
drop temp

*now total from missing agencies
bysort Id2: egen temp = sum(actual_index_total) if total_missing > 0
bysort Id2: egen countymissing = mean(temp)
drop temp

*percent missing
gen propmissing = countymissing/countytotal
replace propmissing = 0 if countytotal == countynonmissing

*create tag
bysort Id2: gen tag = 1 if _n == 1


**
*
*Now Keep Just Never-Missing Agencies and Counties Where Those Account for
* > 90% of crimes
*
**

**Keep 225 counties where never-missing agencies account for more than 90% of crimes
keep if propmissing < .1

**Take a look at total number of crimes
egen temp = sum(actual_index_total)
su temp, f
assert r(mean) == 25325656
drop temp

**Now look at number of crimes omitting those agencies
egen temp = sum(actual_index_total) if total_missing == 0
su temp, f
assert r(mean) == 24922598
drop temp

**This accounts for more than 98 percent of all crimes
disp 24922598/25325656 

**Keep only never-missing agencies
keep if total_missing == 0


**
*
*Collapse
*
**

*Collapse to county month level
collapse (first) agency_name-state_abb agency_type crosswalk_agency_name ///
census_name population_group country_division last_update-agency_count countynonmissing ///
(sum) population_1 officers_killed_by_felony-unfound_index_total, ///
by(fips_state_county_code month firstmonth state county policymonth monthlysanctuary)

*Confirm that collapse worked correctly
bysort fips: egen countytotal = sum(actual_index_total)
assert countynonmissing == countytotal
drop countynonmissing


*-------------------------------------------------------------------------------
* Add in ACS Population Data
*-------------------------------------------------------------------------------

*prep
gen year = 1960 + floor(month/12)
gen fips = fips_s

*merge

merge m:1 fips year using demographics

*drop counties not used
drop if _m == 2

*now all match
assert _m ==3
drop _m


*-------------------------------------------------------------------------------
* Prep variables for analysis
*-------------------------------------------------------------------------------

*Make logs and per capita
foreach var of varlist actual_index_violent actual_index_property actual_index_total {

gen log`var' = ln(`var'+1)
gen `var'percap =(`var'/pop_total) * 100000 * 12
}

*Make clearance rate
gen clearancerate = tot_clr_index_total/actual_index_total

*Make numeric statecounty variable
egen statecounty = group(fips)

*Are there any places with zeros in months? No
bysort statecounty: egen min = min(actual_index_total)
assert min > 0

*-------------------------------------------------------------------------------
* Descriptives
*-------------------------------------------------------------------------------

assert actual_index_violent == actual_rape_total  + actual_robbery_total + actual_assault_aggravated + actual_murder 

*How many total crimes 2013-2015
egen propertysum = sum(actual_index_property) if month >= m(2013m1) & month <= m(2015m12) & monthlysanctuary == 1
su propertysu, detail
assert r(mean) == 4391667 

**
*Table 3
**

latabstat actual_index_violentpercap actual_index_propertypercap ///
actual_index_totalpercap clearancerate, stat(mean median sd max n) tf(table3.tex) f(%9.2f) replace

**California

twoway (lpoly actual_index_totalpercap month if policymonth == m(2013m10) & month > m(2010m1) ///
& month < m(2017m1), yaxis(1) bwidth(1)) ///
(lpoly actual_index_totalpercap month if policymonth > m(2016m10) & month > m(2010m1) ///
& month < m(2017m1), yaxis(1) bwidth(1)), ///
xti("Year") yti("Monthly Annualized Crime Rate (Per 100,000)", axis(1)) ///
xlabel(600(12)672, noticks format(%tmCY)) ///
xline(645, lpattern(dot)) xline(685, lpattern(dot)) ///
legend(order(1  "California" 2 "Nonsanctuary Jurisdictions")) ///
text(3750 645 "CA Sanctuary Law", placement(v) size(small)) 
graph export "S11.pdf", replace



*-------------------------------------------------------------------------------
* Panel Regressions
*-------------------------------------------------------------------------------


***
*Table 4
***

_eststo clear

foreach var of varlist  actual_index_propertypercap actual_index_violentpercap clearancerate{

eststo: areg `var' i.month monthlysanctuary, absorb(statecounty) ///
cluster(state) 

}


esttab using "table4.tex", se mti("Property Crime" "Violent Crime" "Clearance") ///
 ti("Difference-in-Difference Estimates: Effect of Sanctuary on Crime") /// 
coef(monthlysanctuary "Sanctuary" _cons "Intercept") ///
addn("Includes County and Month Fixed Effects; Standard Errors Clustered on State") ///
keep(monthlysanctuary _cons) replace


*-------------------------------------------------------------------------------
* Event Study Regressions
*-------------------------------------------------------------------------------


**
*Make indicator variables
**


gen centeredmonth = month-policymonth

replace centeredmonth = 0 if policymonth == .

*Indicator variables for each centered month variable (only 5 years before and after)
tab centeredmonth if centeredmonth >= -10 & centeredmonth <= 10, gen(indcentered)

* Rename these 
forv j = 1/10{
	loc k = abs(`j' - 11)
	rename indcentered`j' indcenteredN`k'
}
forv j = 11/21{
	loc k = `j' - 11
	rename indcentered`j' indcentered`k'
}

*Now set first indicator to include all previous obs and last to include all subsequent
replace indcenteredN10 = 1 if centeredmonth <= -10
replace indcentered10 = 1 if centeredmonth >= 10 

* Change order of indcenteredN0 (for which y - y* = 0) so that is omitted category
order indcentered0, last

*Replace missing months with zeros 
foreach var of varlist ind* {
replace `var' = 0 if `var' != 1
}


***
*Regressions for Full Period, Excluding CA
***

**Property Crime 

areg actual_index_propertypercap i.month indcenteredN10-indcentered10 if state != "california", absorb(statecounty) ///
cluster(state) 

preserve

* Get coef
regsave  indcenteredN10-indcentered10

* Rename
replace var = "All"

* 95% CI
gen ub = coef + 1.96 * stderr 
gen lb = coef - 1.96 * stderr 

* Generate time variable and other graphing stuff
gen time = _n-11
replace time = time+1 if time > -1

twoway (rspike ub lb time , lw(0.3) lc(gray)) (scatter coef time , msize(small) msymbol(circle) ), ///
	yli(0, lpattern(dash) lc(black)) ///
	xlab(-10(2)10, labsize(small)) ylab(-300(100)200, labsize(small) angle(horizontal)) ///
	title("Property Crime (Excluding California)") ///
	ytitle("") xtitle("") ///
	graphregion( fcolor(white) lcolor(white)) plotregion(fcolor(white) lstyle(none) lcolor(white) ilstyle(none)) ///
	legend(off) xsize(5) ysize(3) name(a, replace)

restore


**Violent Crime 

areg actual_index_violentpercap i.month indcenteredN10-indcentered10 if state != "california", absorb(statecounty) ///
cluster(state) 

preserve

* Get coef
regsave  indcenteredN10-indcentered10

* Rename
replace var = "All"

* 95% CI
gen ub = coef + 1.96 * stderr 
gen lb = coef - 1.96 * stderr 

* Generate time variable and other graphing stuff
gen time = _n-11
replace time = time+1 if time > -1

twoway (rspike ub lb time , lw(0.3) lc(gray)) (scatter coef time , msize(small) msymbol(circle) ), ///
	yli(0, lpattern(dash) lc(black)) ///
	xlab(-10(2)10, labsize(small)) ///
	title("Violent Crime (Excluding California)") ///
	ytitle("") xtitle("") ///
	graphregion( fcolor(white) lcolor(white)) plotregion(fcolor(white) lstyle(none) lcolor(white) ilstyle(none)) ///
	legend(off) xsize(5) ysize(3) name(b, replace)

restore

***
*Regressions for Full Period, Including CA
***

**Property Crime 

areg actual_index_propertypercap i.month indcenteredN10-indcentered10 , absorb(statecounty) ///
cluster(state) 

preserve

* Get coef
regsave  indcenteredN10-indcentered10

* Rename
replace var = "All"

* 95% CI
gen ub = coef + 1.96 * stderr 
gen lb = coef - 1.96 * stderr 

* Generate time variable and other graphing stuff
gen time = _n-11
replace time = time+1 if time > -1

twoway (rspike ub lb time , lw(0.3) lc(gray)) (scatter coef time , msize(small) msymbol(circle) ), ///
	yli(0, lpattern(dash) lc(black)) ///
	xlab(-10(2)10, labsize(small)) ylab(-300(100)200, labsize(small) angle(horizontal)) ///
	title("Property Crime") ///
	ytitle("") xtitle("") ///
	graphregion( fcolor(white) lcolor(white)) plotregion(fcolor(white) lstyle(none) lcolor(white) ilstyle(none)) ///
	legend(off) xsize(5) ysize(3) name(c, replace)

restore


**Violent Crime 

areg actual_index_violentpercap i.month indcenteredN10-indcentered10, absorb(statecounty) ///
cluster(state) 

preserve

* Get coef
regsave  indcenteredN10-indcentered10

* Rename
replace var = "All"

* 95% CI
gen ub = coef + 1.96 * stderr 
gen lb = coef - 1.96 * stderr 

* Generate time variable and other graphing stuff
gen time = _n-11
replace time = time+1 if time > -1

twoway (rspike ub lb time , lw(0.3) lc(gray)) (scatter coef time , msize(small) msymbol(circle) ), ///
	yli(0, lpattern(dash) lc(black)) ///
	xlab(-10(2)10, labsize(small)) ///
	title("Violent Crime") ///
	ytitle("") xtitle("") ///
	graphregion( fcolor(white) lcolor(white)) plotregion(fcolor(white) lstyle(none) lcolor(white) ilstyle(none)) ///
	legend(off) xsize(5) ysize(3) name(d, replace)

restore

***
*
*Figure S12
*
***


graph combine c d a b, xsize(3) ysize(2) ///
l1title("Change in Crime Rate Relative to Date of Policy", size(med)) ///
b1title("Months Since Sanctuary Policy")

graph export "S12.pdf", replace





*-------------------------------------------------------------------------------
* Event Study Regressions, Balanced Panels
*-------------------------------------------------------------------------------












***
*Regressions for Early Balanced Panel
***

preserve

**
*Keep Early Panel Only


keep if firstmonth <= m(2010m12)
keep if month >= m(2010m12) 

*Ten month window on either side of treatment
gen ten = 1 if centeredmonth == 10 | centeredmonth == -10

egen ten_window = sum(ten), by(statecounty)

keep if ten_window == 2

count
assert r(N) == 3172

**Number of observations
codebook statecounty
count

**Property Crime 

areg actual_index_propertypercap i.month indcenteredN10-indcentered10 , absorb(statecounty) ///
cluster(state) 

* Get coef
regsave  indcenteredN10-indcentered10

* Rename
replace var = "All"

* 95% CI
gen ub = coef + 1.96 * stderr 
gen lb = coef - 1.96 * stderr 

* Generate time variable and other graphing stuff
gen time = _n-11
replace time = time+1 if time > -1

twoway (rspike ub lb time , lw(0.3) lc(gray)) (scatter coef time , msize(small) msymbol(circle) ), ///
	yli(0, lpattern(dash) lc(black)) ///
	xlab(-10(2)10, labsize(small)) ///
	title("Property Crime (Early Panel)") ///
	ytitle("") xtitle("") ///
	graphregion( fcolor(white) lcolor(white)) plotregion(fcolor(white) lstyle(none) lcolor(white) ilstyle(none)) ///
	legend(off) xsize(5) ysize(3) name(a, replace)

restore

preserve

**
*Keep Early Panel Only

keep if firstmonth <= m(2010m12)
keep if month >= m(2010m12) 

*Ten month window on either side of treatment
gen ten = 1 if centeredmonth == 10 | centeredmonth == -10

egen ten_window = sum(ten), by(statecounty)

keep if ten_window == 2


**Violent Crime 

areg actual_index_violentpercap i.month indcenteredN10-indcentered10, absorb(statecounty) ///
cluster(state) 

* Get coef
regsave  indcenteredN10-indcentered10

* Rename
replace var = "All"

* 95% CI
gen ub = coef + 1.96 * stderr 
gen lb = coef - 1.96 * stderr 

* Generate time variable and other graphing stuff
gen time = _n-11
replace time = time+1 if time > -1

twoway (rspike ub lb time , lw(0.3) lc(gray)) (scatter coef time , msize(small) msymbol(circle) ), ///
	yli(0, lpattern(dash) lc(black)) ///
	xlab(-10(2)10, labsize(small)) ///
	title("Violent Crime (Early Panel)") ///
	ytitle("") xtitle("") ///
	graphregion( fcolor(white) lcolor(white)) plotregion(fcolor(white) lstyle(none) lcolor(white) ilstyle(none)) ///
	legend(off) xsize(5) ysize(3) name(b, replace)

restore



***
*Regressions for Late Balanced Panel
***

preserve

**
*Keep Late Panel Only

keep if firstmonth <= m(2012m8)
keep if month >= m(2012m8) 

*Ten month window on either side of treatment
gen ten = 1 if centeredmonth == 10 | centeredmonth == -10

egen ten_window = sum(ten), by(statecounty)

keep if ten_window == 2


**Number of observations
codebook statecounty
count

**Property Crime 

areg actual_index_propertypercap i.month indcenteredN10-indcentered10 , absorb(statecounty) ///
cluster(state) 

* Get coef
regsave  indcenteredN10-indcentered10

* Rename
replace var = "All"

* 95% CI
gen ub = coef + 1.96 * stderr 
gen lb = coef - 1.96 * stderr 

* Generate time variable and other graphing stuff
gen time = _n-11
replace time = time+1 if time > -1

twoway (rspike ub lb time , lw(0.3) lc(gray)) (scatter coef time , msize(small) msymbol(circle) ), ///
	yli(0, lpattern(dash) lc(black)) ///
	xlab(-10(2)10, labsize(small)) ///
	title("Property Crime (Late Panel)") ///
	ytitle("") xtitle("") ///
	graphregion( fcolor(white) lcolor(white)) plotregion(fcolor(white) lstyle(none) lcolor(white) ilstyle(none)) ///
	legend(off) xsize(5) ysize(3) name(c, replace)

restore

preserve

**
*Keep Late Panel Only

keep if firstmonth <= m(2012m8)
keep if month >= m(2012m8) 

*Ten month window on either side of treatment
gen ten = 1 if centeredmonth == 10 | centeredmonth == -10

egen ten_window = sum(ten), by(statecounty)

keep if ten_window == 2


**Violent Crime 

areg actual_index_violentpercap i.month indcenteredN10-indcentered10, absorb(statecounty) ///
cluster(state) 

* Get coef
regsave  indcenteredN10-indcentered10

* Rename
replace var = "All"

* 95% CI
gen ub = coef + 1.96 * stderr 
gen lb = coef - 1.96 * stderr 

* Generate time variable and other graphing stuff
gen time = _n-11
replace time = time+1 if time > -1

twoway (rspike ub lb time , lw(0.3) lc(gray)) (scatter coef time , msize(small) msymbol(circle) ), ///
	yli(0, lpattern(dash) lc(black)) ///
	xlab(-10(2)10, labsize(small)) ///
	title("Violent Crime (Late Panel)") ///
	ytitle("") xtitle("") ///
	graphregion( fcolor(white) lcolor(white)) plotregion(fcolor(white) lstyle(none) lcolor(white) ilstyle(none)) ///
	legend(off) xsize(5) ysize(3) name(d, replace)

restore




graph combine a b c d, xsize(3) ysize(2) ///
l1title("Change in Crime Rate Relative to Date of Policy", size(med)) ///
b1title("Months Since Sanctuary Policy")

graph export "S13.pdf", replace


















